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Abstract — In the last two decades the village of Casalabate was affected by a number of sinkhole events, the last episode 
occurring on 7 January 2012. Due to some geological and geotechnical uncertainties, Casalabate became an intriguing case 
study of sinkhole hazard in urban areas and three geological-geophysical models were recently proposed by three different 
studies. In this paper new geophysical data and some epistemological consideration were done. The results are a 
contribution to an efficient predictive model that could be well define both the geological setting and sinkhole hazard. 
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I. Introduction 

Casalabate is a coastal small town frequently affected by sinkhole events. The first noticed event occurred on 1993 and 
caused several destructions of buildings. In the same year, a sinking of a small portion of emerged beach happened (Delle 
Rose and Federico, 2002). Four years later a collapse affected the submerged portion of beach and allowed the formations of 
a submarine spring, since then the bathing was forbidden for some weeks. Two episodes of sinkhole on the emerged beach 
occurred during 2000 and 2004, along a strike no longer than 30 meters. Each of these caused on surface short-lived elliptical 
depression quickly filled by sand. All the aforementioned episodes occurred inside a narrow area (about 250 m long and 100 
m wide) close to the coast and oriented according to a north-north-east - south-south-west tectonic structures system (Delle 
Rose and Leucci, 2010). Further episodes occurred on March 2010, November 2010 and March 2011. The last noticed 
sinkhole happened on 7 January 2012. Usually the sinkholes occurred during or after rains or sea storms and each episode 
lasts some minutes at least. Due to geological and geotechnical uncertainties, Casalabate became an intriguing case study of 
sinkhole hazard in urban areas and three geological-geophysical models were recently proposed to explain the ground- 
instability and to asses the related hazard (Delle Rose and Leucci, 2010; Calo et al., 2011; Margiotta et al., 2012). 3D 
geological model and new geophysical data here shown can contribute to an efficient predictive model taking into account 
the hydrogeology context. Moreover, some epistemological considerations had performed to improve the interdisciplinary 
exchange between geologists and geophysicists. 

II. Geological setting 

Casalabate is located on the Adriatic side of the Salento Peninsula (Southern Italy) on a coastal plain elevated a few meters 
above sea level (Fig. 1). The Salento Peninsula is the emerged southeastern portion of the Adriatic Carbonatic Plate (Auct.) 
which is formed by Jurassic-Cretaceous limestones and dolostones covered by Tertiary and Quaternary carbonates and marls. 
The bedrock of Casalabate is build up by Quaternary carbonate deposits, some tens of meters thick, lie on Cretaceous 
limestone (Rossi, 1969). Nevertheless, geological and geotechnical data suitable for civil engineering lacked until few years 
ago. As a consequence, to face the sinkhole hazard the municipal authority performed an extensive bedrock survey by means 
of boreholes and geophysical measures. The results of such a survey were used in Delle Rose and Leucci (2010) and in Calo 
et al. (2011). Casalabate bedrock resulted built up by “two overlap units”. The lower of which is a moderately cemented 
massive calcarenite made up mainly of marine bioclasts and medium-large in grain size. The upper units is formed by sets of 
variableness cemented thinly laminated and cross-stratified calcarenite [ancient dunes] whose grains were terrestrial or 
marine supplied. The latter is 5-10 m thick, whereas the base of the former are locally unknown but, on the base of the setting 
surrounding Casalabate, the thickness of the massive calcarenite doesn’t excess 20-25 m. Both the calcarenite units show a 
nearly omogeneous low-angle beds dip, which seems to be due to primary clinostratification as well as tectonic eastward 
tilting” (Delle Rose and Leucci, 2010). Four tectonic fracture systems (medium strikes N 5°E, N55°E, N100°E and N145°E) 
probably affect both the stratigraphic units. The shape and orientation of the landforms are clearly due to the spatial 
arrangement of the tectonic fractures. This structural control is particularly evident as it concerns the shape of the depressions 
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filled by clayey-silt marsh deposits. The calcarenite units were interested by various periods of dissolution karst processes to 
begin from the Middle Pleistocene. Usually such processes are named as “parakarstic” (Auct.) to emphasize the small 
dimension of the morphological forms respect those charactering the Jurassic -Cretaceous limestones and dolostones. 



Fig. 1 - Surveyed area (photo by Google Earth) 

A shallow coastal aquifer flows out along the Casalabate coast. The groundwater circulation is influenced by the tectonic 
fracture systems and is manifested by some submarine springs. Underwater investigations of the sinking of a portion of 
submerged beach in 1997 revealed the presence of a karst spring (Delle Rose and Federico, 2002) probably joined with the 
cave that induced the 1993 sinkholes. Such springs are not perennial and can occur suddenly. The last case of new spring 
formation occurred below an house on the beginning of September 2013 and caused a new warning to the residents. 

III. 3D Geological model 

Using the software Groundwater Modeling System (GMS), a 3D geological model and a 3D groundwater model are in 
progress. The scope of building the models is to simplify the field problem and organize the field data so that the system can 
be analyzed more readily (Anderson and Woessner, 2002). The conceptualization of the models includes synthesis and 
framing up of data pertaining to geology, hydrogeology, hydrology, and meteorology. The lithologic data were collected 
from 23 core lithologs data (kindly provided by the technical office of Lecce municipality). Such a borehole survey was 
performed on 2008 and was related to the project “actions to protect the soil in the hydrogeological risk areas of the 
municipality of Lecce” funded by the Ministry of Environment and Protection of Land. 

The geological model was created using a block-centric finite -difference grid. The gridding was performed through the eight 
nearest-neighbor methodology with 3-D interpolation by average minimum distance. The solid model was developed to a 
depth of about 10 m below the ground level. The sensitivity of the model was tested by varying the horizontal and vertical 
spacing of the nodes, and an optimized final model, least sensitive to changes in spatial resolution (i.e., smaller grid sizes), 
was built. The resolution of the final model was 100mxl00mx2m. The resulting discretization consisted of 
250nodesx250nodesxl00nodes, obtaining 6,250,000 solid model nodes, each with a voxel volume of 20,000 m 3 . It should be 
noticed that the results of the modeling are not free from uncertainties, which could be reduced by increasing data points, but 
they illustrate one of the most probable scenarios. The smallest scale of variation that the model is able to depict is equal to 
the resolution of one voxel. The 3D geological model (Fig. 2a) with the cross-sections in Fig. 2b provides a depictions of the 
bedrock of the studied area. Such model was vertically discretized in 8 layers (paving; artificial landfills; sands of beach; 
organic soils; residual soils; marsh deposits; cemented thinly laminated and cross-stratified calcarenites; moderately 
cemented massive calcarenites). The layers were allowed to have seepage from the top and leakage through the base, making 
them hydraulically connected. The top-most layer and the second layer were defined as unconfined. The rest of the layers 
were defined as confined, as the water table was not expected to fall >10 m. Such a model has reached a first level of 
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advancement and needs several new data to be fully efficient, to start from the boundary conditions as well as the real 
physical characters of the higher layers. 
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Fig. 2 - 3D Geological model (a) and cross sections ol Casalabate urban area (b) 

Given the heterogeneity and complexity associated with the study area (e.g. local-scale variations of aquifer properties), 
some generalizations, simplifications and assumptions will be made to construct the groundwater flow model. The modeling 
area is limited in the east by the Adriatic sea, and it accounts for about 1,5 km in length and 0,5 km in width. The 3-D 
constant-density groundwater flow will be simulated by a block-centric, finite -difference grid model, using MODFLOW 
code inside GSM software. Furthermore MODAEM code, which allows the analysis to the analytical elements (polygons, 
lines and points), will be used. 

IV. New geophysical measurements 

In our previous work, we have proposed a method to integrate geophysical survey and geological knowledge to characterize 

Page | 73 
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the sinkhole hazard at Casalabate. As regards geophysics work, the performed 3D Ground Penetrating Radar (GPR) and 3D 
Electrical Resistivity Tomography (ERT) allowed us to reveal a complex hypogean cavity system inside the urban area 
affected by sinkhole since the 1993 (Delle Rose and Leucci, 2010). Here we propose ERT and GPR survey carried out in an 
area (Fig. 3) next to the 1993 sinkhole collapsed that caused the collapsing of three buildings and irreparable damage to 
various others, which were demolished later. 



> 2D ERT profile 
I 1 3D ERT surveyed area 

I I 3D GPR surveyed area 

r - ' Collapsed area 

v / 

FIG. 3 - ERT AND GPR SURVEYED AREAS 

For ERT measurements a 48-channel Syscal-Pro Resistivity-meter in a multi electrode configuration using a Wenner— 
Schlumberger array (Loke 2009) was used. This is a hybrid between the Wenner and Schlumberger arrays (Pazdirek and 
Blaha 1996) arising out of more recent work with electrical imaging surveys. This array is moderately sensitive to both 
horizontal and vertical structures. In areas where both types of structures are expected, this array might be a good 
compromise between the Wenner and the dipole -dipole array (Loke 2009). 

GPR measurements were performed using a georadar Ris Hi Mod with a 600 MHz (centre frequency) antenna manufactured 
by Ingegneria dei Sistemi (IDS). 

V. 2D ERT Measurements 

The 48 electrodes were located on the road surface and spaced 5 m. The line were oriented predominantly in a W-E direction. 
The measurements were stacked at least three times for each point along the profile to enhance data quality. Standard 
deviation was used for evaluating the statistical dispersion of recorded resistivity data from their arithmetic mean. The 
chosen electrode array and spacing made it possible to investigate about 30 m of subsoil with a good resolution. Resistivity 
data were inverted with the TomoLab software (Geostudi Astier, 2010). The 2D resistivity model is shown in figure 4. It 
shows the resistivity distribution of the upper 30 m, in which it is possible to note the following: (i) a layer of relatively low 
resistivity (about 50 ohm m), from the surface to about 2 m in depth (labeled A); (ii) a layer (labeled B) in which resistivity 
values range between 100 and 150 ohm m; (iii) the area of highest resistivity (p = 200 ohm m), inside the layer B; iv) a 
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deeper layer (labeled C) in which resistivity values range between 100 and 130 ohm m. The model established by means of 
2D resistivity imaging profiles thus allows three different zones to be detected and analysed. The first one (A) may 
correspond to the anthropogenic layer (paving; artificial landfills) in which a groundwater with salt water intrusion is present. 
The second one (B) corresponds to lamp strata of the laminated and cross-stratified calcarenites and the third one (C) to the 
moderately cemented massive calcarenites. 



Fig. 4- 2D ERT PROFILE: RESISTIVITY distribution in the subsoil 

The area of highest resistivity inside the layer B could be related to karstic forms. The location of such a forms well fits a 
complex marsh depression (probably a group of coalescent dolines or sinkholes) buried below the anthropogenic layer that 
was detected by means of photo analysis of aerial photo taken before urbanization (Delle Rose and Fiorito, 2000, Fig. 1 of p. 
76). 

VI. 3D ERT Measurements 

Twenty ERT lines were acquired at the site over an area of 45x45m 2 . The electrode spacing was of 1 m. The ERT line 
spacing was irregular due to the presence of several buildings. A total of 25500 data points were measured. To avoid damage 
to the road surface, electrocardiogram electrodes were used. As demonstrated by Leucci (2006, 2007), the electrical 
resistance contact (ERC) value was of about 2.0 kQ. As affirmed by Zonge and Hughes (1985), the results were unaffected 
by this ERC value. The voltage value was set at max voltage, while the injected current values were automatically adjusted 
by the recording system. Then the measurements were stacked at least three times for each point along the profile to enhance 
data quality. Standard deviation was used for evaluating the statistical dispersion of recorded resistivity data from their 
arithmetic mean. 

The acquired data were combined to produce three dimensional volumes of the resistivity distribution in the subsoil. The 
ErtLab software manufactured by Multi-Phase Technologies, LLC, was used to automatically invert the resistivity acquired 
data and to yield a 3D rho model. Its numerical core based on tetrahedral FEM and inversion was performed with a robust 
inversion (data variance iterative reweighting). (Geostudi Astier, 2010). The results of the inversion of the electrical data set 
are shown in Fig 5a and b. 

It shows the resistivity distribution in which it is possible to note the following: (i) a layer of relatively low resistivity (about 
40-50 ohm m), at depth from the surface ranging from 2m to 3m in depth (labelled A); (ii) a layer (labelled B) in which 
resistivity values range between 100 and 110 ohm m; (iii) a deeper layer (labeled C) in which resistivity values range 
between 70 and 80 ohm m; iv) the area of highest resistivity (p = 200 ohm m) labeled D. The model established by means of 
3D resistivity imaging profiles thus allows three different zones to be detected and analysed. The three zones are similar to 
those described in the 2D ERT profile analysis. Also here the area of highest resistivity could be related to karstic forms. 
Particularly the zone labeled “V” could be related to a developed void zone. 

The 3D images of resistivity can easily be visualized by 3D contouring of iso -resistivity volumes (Fig. 6). In this 
representation, the transparency function is defined by two threshold values of the rho, rhol and rho2 (rhol < rho2). In the 
intervals rho < rhol and rho > rho2, data are rendered as transparent, therefore only the data in the interval rhol < rho < rho2 
are visualized. The threshold calibration is a very delicate task. In fact, by lowering the threshold value, not only the visibility 
of the main anomaly is raised, but also that of the smaller objects and noise increases. In Fig. 6, the Ip data set is displayed 
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with iso-IP volumes using a threshold value ranging from 80 to 200Qm. The continuous high rho anomaly labeled V is more 
visible. This kind of visualization allows the anomalies resistivity, already shown in Fig. 5, to be emphasized. It is clear the 
spatial position of karstik form. 



FIG. 5- 3D ERT VOLUME ANALYSIS 
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Fig. 6: 3D Contouring of iso-resistivity volumes 
VII. 3D GPR Measurements 

A reconnaissance survey was conducted in a continuous mode in an area of 55 m x 12 m. Due to the presence of several 
buildings, GPR profiles were acquired in a grid with variable dimensions ranging from 0.5m to lm. The quality of the raw 
data did not require advanced processing techniques. However, appropriate processing was performed for easier 
interpretation. The processing steps were a background removal filter process (Conyers and Goodman 1997) and 2D velocity 
migration (Yilmaz 1987). In order to transform the time axis in depth axis the electromagnetic (EM) wave velocity was 
esteemated in the surveyed area. EM wave velocity can be estimated from GPR data in several ways; the conventional 
method involves common depth-point (CDP) and wide-angle reflection and refraction (WARR) data sets. Both methods 
require two antennas in separate units and relatively long acquisition times. In the first case, both antennas are 
simultaneously moved apart on either side of the mid-point of 

the profile. In the second case, one antenna remains stationary while the other is moved along the profile direction. The EM 
wave velocity can be more quickly and easily determined from the reflection profiles acquired in a continuous mode, using 
the characteristic hyperbolic shape of reflection from a 

point source (diffraction hyperbola). This is a very common method for the EM velocity estimation and is based on the 
phenomenon that a small object (the dimensions of the object are smaller than the wavelength of the EM wave introduced in 
the ground) reflects EM waves in almost every direction. 

This method is supported by most radar software packages and should be applied if the diffraction hyperbola occurs at least 
in 20 radar scans (Fruhwirth and Schmoller 1996). In the data acquired several hyperbolic reflections (Fig. 7), which allow 
an accurate velocity analysis, are present. The velocity analysis was performed using the GPR-Slice software (Goodman, 
2013), on the radar data set, and the application of this method points out an average EM wave propagation velocity of 0.07 
m ns '. Figure 7 shows one of the profiles acquired after the above described processing. A general characteristic of the 
surveyed area seems to allow a good penetration of electromagnetic (EM) energy (about 60 ns corresponding to a depth of 
about 2.1 m if the mean velocity value of 0.07 m ns-1 is used). This is essentially due to the physical characteristics of the 
subsurface, made up of material which poorly dissipates EM energy. This is also the case in all the other profiles acquired in 


an area. 
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In the processed radar sections, several reflection events are clearly identifiable. At times ranging from about 20 ns to about 
40 ns (0.7 m to 1 .4 m in depth) it is possible to identify some reflection events likely due to the probable of karstic forms . 



Fig. 7: GPR processed profile 

Since GPR data were acquired in a grid of closely spaced profiles, reflection amplitude can be accurately placed in three 
dimensions. A way to obtain visually useful maps for understanding the plan distribution of reflection amplitudes within 
specific time intervals is the creation of horizontal time slices. There are maps on which the reflection amplitudes have been 
projected at specified time (or depth), with a selected time interval (Conyers, 2006). In a graphic method developed by 
Goodman et al., (2006), named “overlay analysis”, the strongest and weakest reflectors at the depth of each slice are assigned 
specific colors. This technique allows the linkage of structures buried at different depths. This represents an improvement in 
imaging because subtle features that are indistinguishable on radargrams can be seen and interpreted in an easier manner. In 
the present work the time slice technique has been used to display the amplitude variations within consecutive time windows 
of width At = 5 ns. The time slices (Fig. 8) show the normalized amplitude using a range that define the blue color as the zero 
level and red color as 1 level. Fig. 8 shows the most significant slices. In the slices ranging from 0.7m to 1.36m depth, 
relatively high amplitude anomalies are clearly visible in the dashed dark box. These correspond to the karstic forms areas. 
Moreover the highest amplitudes were rendered into an isosurface (Conyers, 2004, Conyers, 2012, Conyers et al., 2013). 3D 
amplitude isosurface rendering displays amplitudes of equal value in the GPR study volume. Shading is usually used to 
illuminate these surfaces, giving the appearance of real archaeological structures. In this case the threshold calibration is a 
very delicate task in order to obtain useful results. 
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This visualization technique better evidence the anomalies found in the surveyed area. In this case is clearly visible the 3D 
development of the anomalies related to the karstic forms (Fig. 9). 



Fig. 9: GPR data 3D visualization by the EM wave amplitude iso-surface 


VIII. How GEOLOGICAL UNCERTAINTIES WORSEN THE HAZARD ASSESSMENT 

How previous said, three geological-geophysical models were proposed to explain the instability of the underground of 
Casalabate (Delle Rose and Leucci, 2010; Calo et al., 2011; Margiotta et al., 2012). The first of these considers a shallow 
karst cavities system in carbonate rocks filled by different kind of deposits (marsh deposits; residual soils; organic soils; 
sands of beach) suitable to be removed by water flows. In such a context, sinkholes could be due to both the break of 
cavities’ vault (i.e. collapse sinkhole, cf. Waltham et al. 2005) and the down-washing of the filling deposits inside the cavities 
system (i.e. suffosion sinkhole, cf. Waltham et al. 2005). This latter phenomena “is particularly probable taking into account 
both the sea level changes and the variations of the piezometric level of the groundwater inside the hypogean system, that can 
lead high energy water flows that are able to remove and wash away the sediment filling” (Delle Rose and Leucci, 2010, p. 
153). According to Calo et al. (2011) the sinkhole “phenomena cannot be due to the collapse of the underground cavities 
vault, but because of the incidence and combination of the following factors: poor, and in some cases, very poor geotechnical 
characteristics of foundation grounds; the lack of a preliminary geotechnical characterization of buildings leading to the 
choice of improper practice in laying down foundations; the letting of rainwater into the subsoil which caused the leaching 
and the removal of the fine fraction of the foundation grounds”. Finally, Margiotta et al. (2012) claim the occurrence of 
collapse sinkholes inside the marsh deposits and suffosion sinkholes moving down sandy beach deposits. Such these Authors 
emphasize how their “data” with “particular regard to the stratigraphic al succession, are in strong disagreement” with the 
model of Delle Rose and Leucci (2010) whereas in the same time they are suitable to “mapping the susceptibility to 
sinkhole” at Casalabate. So, apparently these three models are incompatible each with another one. Really the divergences 
are mainly due to the different methods to achieve the scientific knowledge or, in other words, they concern epistemological 
aspects of the research. 

Both the models proposed by Delle Rose and Leucci (2010) and Calo et al. (2011) used a stratigraphy reconstruction based 
the boreholes drilled in the aforementioned project “actions to protect the soil in the hydrogeological risk areas”. However, 
“outcrops available to geological observations are scarce in the Casalabate area, so a borehole survey and sampling had been 
the key tools to reconstruct the main features of the local geological substratum” (Delle Rose and Leucci, 2010). The subsoil 
detection performed let us to recognize in the first ten meters of the bedrock two overlap units, the lower of which is a 
moderately cemented massive calcarenite and the upper unit is formed by sets of variably cemented thinly laminated and 
cross-stratified calcarenite. As regard the differences between the model proposed by Delle Rose and Leucci (2010) and Calo 
et al. (2011), these last Authors stressed how “the integrated analysis of information obtained from logs and georadar 
investigation evidences that [...] no cavities and/or channels referable to parakarstic dissolution phenomena have been 
detected”. We point out that only 2 out of 23 boreholes performed on 2008 were performed within the area really affected by 
sinkholes (the “collapsed area” evidenced in Fig. 3a). Moreover, one of the such two boreholes crossed red residual soils that 
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are a typical filling of dolines. So, in our opinion, the conclusion of Calo et al. (2011) can be relevant at the entire urban scale 
but is inconsistent with the scale of the dangerous area. As before said such an area where sinkholes occurred is narrow 
(about 250 m long and 100 m wide) and close to the coast. However we have to observe that Calo et al. (2011) did not show 
the profiles of their georadar investigation and therefore don’t allow a critical evaluation of their measures and data. 

The divergences between the models of Delle Rose and Leucci (2010) and Margiotta et al. (2012) might be explained from 
the different results of the respective stratigraphy and geophysical surveys. It must be observed that Margiotta et al. (2012) 
ignored both the aforementioned borehole survey performed on 2008 and the paper of Calo et al. (2011). However they 
misunderstand the model of Delle Rose and Leucci (2010) especially when they assert that such a model ‘‘is based on a 
stratigraphy consisting exclusively of carbonate deposits, without taking into the due account the crucial presence of cover 
deposits above the calcarenite bedrock, and primarily that of the organic -rich clays [i.e. marsh deposits]” (Margiotta et al., 
2012, p. 667). 

Margiotta et al. (2012) described the marsh deposits with a thickness of up to 11 m and overlying (without interruptions) the 
calcarenites along the coast for an outcropping area ranging from 100 to 250 meters. This stratigraphy setting does not fit the 
boreholes data of the aforementioned project funded by the Ministry of Environment and Protection of Land: in fact only 6 
perforations out of 23 crossed palustrine deposits of thicknesses ranging from 0.3 to 1.4 m. Moreover, the whole 
stratigraphical model proposed in Margiotta et al. (2012) appears based on subjective considerations. The Authors affirmed 
that their “research involved [. . .] the critical revision of more than 40 well-core stratigraphies”, without describing the nature 
of such a “critical revision”. Their subsoil model seems mainly based on interpretation of geophysical data. Consequently the 
controversy between the models embraced also methodological aspects regarding data processing and interpretation. We 
have to observe that for geophysical data Margiotta et al. (2012) associate individual radar signal reflections to specific 
geological boundaries and even associate a radar signal reflection to a prediction of subsidence. However, it is common 
practice to perform the radar profiles in correspondence of stratigraphical columns or trenches to be able to associate specific 
radar reflections to geological stratigraphy (Conyers and Goodman, 1997; Conyers, 2004). A closer look at the geophysical 
data of Margiotta et al. (2012) is really difficult because the figures were not possible to read the metric scales, time scales, 
and resistivity scales. Their radar sections are surely affected by ringing (Conyers and Goodman, 1997; Conyers, 2004) 
probably related to low penetration of radar signal. Seems that the radar signal do not propagate below 20ns. This is due to 
the high conductive materials presents in the subsoil as confirmed by ERT results presented by the same Authors. In fact in 
the surveyed area electrical resistivity seems to ranging from 0.5 to 20 ohm m, index of highly conductive subsurface. 
Therefore the anomalous zones below 20 ns could probably be an artifact due to the application of the background removal 
filter. For ERT profiles Margiotta et al. (2012) associate resistivity values to specific geological formations. For example 
0.5-2 ohm m was related to clay sediments without any reference to a nearby stratigraphycal core or a bibliography. In fact 
Reynolds, (2011) on page 291, presented a table in which he correlates various types of geological materials with a range of 
values of electrical resistivity (for clay materials resistivity value ranges from 1 to 150 ohm m). Also in Loke (2011, page 6) 
for clay materials, resistivity value ranges from 1 to 100 ohm m, while for the sea water resistivity value is about 0.2 
ohm □ m (due to the relatively high salt content). Given the above, the stratigraphical features of the marsh deposits 
described by Margiotta et al. (2012) appear inconsistent and the result of a wrong interpretation of geophysical data. 
Consequently, their model of sinkhole susceptibility is questionable. 

IX. Conclusion 

The case of Casalabate is a good example of how geological uncertainties can worsen the hazard assessment. In 
consideration of the complexity of the problem and the concerns of the collectivity, an efficient predictive model of the 
hazard should be defined. In order to select the most appropriate geological-geophysical model, confrontations among all the 
researchers working on this issue would be desirable. Especially the knowledge exchange between geologists and 
geophysicists must be improved. As revealed the ERT and GPR surveys exposed in this paper, dangerous cavities could be 
placed also out of the area repeatedly interested by sinkholes since 1993. Consequently an efficient predictive model taking 
into account the hydrogeology and hydraulic context, especially the occurrence of rains and sea storms, must be produced. 
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